package test.demo;

import tech.waterism.modelbase.PredictBase;
import tech.waterism.modelbase.PredictResult;
import com.alibaba.fastjson.JSON;

import java.util.*;

/**
 * 河北雨洪模型 河道演算
 * @author 初京刚
 * @date 2022-4-17
 * @version 1.0
 */

public class HBHDMODEL extends PredictBase {

    /**
     * 河段长（km）
     */
    private double L;

    /**
     * 河宽（m）
     */
    private double W;

    /**
     * 起始下渗率，mm/h
     */
    private double Hfi;

    /**
     * 稳定下渗率，mm/h
     */
    private double Hfc;

    /**
     * 表示土壤物理性质的下渗曲线指数
     */
    private double Hk;

    /**
     * 下渗量计算的允许误差
     */
    private double err_limit_f;

    /**
     * 1km河长的汇流时间（分钟）
     */
    private double Ah;

    /**
     * 河道断面形状参数
     */
    private double bh;

    /**
     * 河段上下断面或流域平均最大雨强发生时间至洪峰出现的传播时间（时段）
     */
    private int Th;

    /**
     * 汇流计算的允许误差
     */
    private double err_limit_q;

    /**
     * 降水时间
     */
    private Date[] floodTm;

    /**
     * 降水量
     */
    private double[] floodDrp;

    /**
     * 降水序列长度
     */
    private Integer floodRainRange;

    /**
     * 降水时段长
     */
    private double timeInterval;

    /**
     * 上游断面流量过程
     */
    private double[] IA;

    /**
     * 时段出流状态
     */
    private double[] Q0;

    /**
     * 时段末出流
     */
    private double[] Q;


    public HBHDMODEL(Map<String, String> data) {
        super();
        L = Double.parseDouble(data.get("L"));
        W = Double.parseDouble(data.get("W"));
        Hfi = Double.parseDouble(data.get("HFI"));
        Hfc = Double.parseDouble(data.get("HFC"));
        Hk = Double.parseDouble(data.get("HK"));
        Ah = Double.parseDouble(data.get("AH"));
        bh = Double.parseDouble(data.get("BH"));
        Th = Integer.parseInt(data.get("TH"));
        floodTm = JSON.parseObject(data.get("dt"), Date[].class);   // input[i].data.dt
        floodDrp = JSON.parseObject(data.get("rain"), double[].class);// input[i].data.rain
        IA = JSON.parseObject(data.get("q"), double[].class);      // input[i].data.q
        timeInterval = Double.parseDouble(data.get("clen"));        // index.clen
        floodRainRange = floodDrp.length;
        Q0 = new double[IA.length + 1];
        Q0[0] = Double.parseDouble(data.get("STARTQ"));
        Q = new double[IA.length];

        err_limit_q = 0.01;
    }

    /**
     * 模型计算
     */
    public Map<String, double[]> HBHDmodel() {

        //汇流计算
        HDSS();

        Map<String, double[]> simulationResult = new HashMap<>();
        simulationResult.put("rain", floodDrp);
        simulationResult.put("QSim", Q);

        return simulationResult;
    }

    /**
     * 有沿程损失的河道汇流计算
     *
     * @return 出流过程
     */
    private void HDSS() {
        double MI = 0.0; //时段累积入流，包括上游入流和河面降雨 m3
        double MF = 0.0; //时段累积下渗 m3
        double IW = 0.0; //时段累计入流与下渗之差
        double IW0 = 0.0; //上一时段累计入流与下渗之差
        double Ia = 0.0; //时段上游平均入流 m3/s
        double Fh = 0.0; //小时下渗量 m3/h
        double fh = 0.0; //下渗率 mm/h
        double Ph = 0.0; //时段河面平均降雨 m3/h
        double IP = 0.0; //时段入流 m3

        double Ut = 3600 * timeInterval; //时间单位转换

        for(int i = 0; i < floodRainRange; i++) {
            MI = 0.0;
            MF = 0.0;
            if(i < Th) {
                Ia = 0.0;
            } else {
                Ia = IA[i - Th];
            }
            Ph = floodDrp[i] * W * L / timeInterval;
            for(int t = 0; t < (int)Math.round(timeInterval); t++) {
                fh = (Hfi - Hfc) * Math.exp((-1) * Hk * t) + Hfc;
                Fh = fh * W * L;
                MI += (Ia * 3600 + Ph);
                MF += Fh;
            }
            IW = MI - MF;
            if(IW > 0.0){
                IP = (IW - IW0)/Ut;
                IW0 = IW;
                HDHL(i, IP);
            } else {
                IW0 = 0.0;
                Q[i] = 0.0;
            }
            Q0[i + 1] = Q[i];
        }
    }

    /**
     * 河道汇流计算
     *
     * @param i 时段
     * @param IAH 入流
     * @return 出流过程
     */
    private void HDHL(int i, double IAH) {
        double Beta = 0.9;
        double AAH = Ah * Math.pow(L, Beta) / 60 / timeInterval;
        double SoH0;//时段初河道蓄水量 m3
        double QoH0;//时段初河道出流 m3/s
        double QMA;//时段末流量迭代求解已知部分
        double QTT;//时段末流量迭代求解中间判断
        double QTT2;//时段末流量迭代求解主体

        double Ut = 3600 * timeInterval; //时间单位转换

        if(i == 0) {
            QoH0 = Q0[0];
        } else {
            QoH0 = Q[i - 1];
        }

        if(IAH > 0.0) {
            if(bh == 0.5) {
                Q[i] = Math.pow(((-2) * AAH + Math.pow(4 * AAH * AAH - QoH0 * Ut * Ut
                        + 4 * AAH * Math.pow(QoH0, 1 / 2) * Ut + 2 * IAH * Ut * Ut, 1 / 2)) / Ut, 2);
            } else if (bh == 1.0) {
                if (QoH0 > 0.0) {
                    SoH0 = AAH * (Math.exp(QoH0) - Math.exp(0.01));
                } else {
                    SoH0 = 0.0;
                }
                QMA = IAH - 0.5 * QoH0 + SoH0 + AAH * Math.exp(0.01);
                QTT = QoH0;
                QTT2 = Math.exp((QMA - 0.5 * QTT) * Ut / AAH);
                while(Math.abs(QTT2 - QTT) > err_limit_q) {
                    QTT = QTT2;
                    QTT2 = Math.exp((QMA - 0.5 * QTT) * Ut / AAH);
                }
                Q[i] = QTT2;
            } else if(bh < 1.0) {
                QMA = IAH - 0.5 * QoH0 + (AAH / ((1 - bh) * Ut)) * Math.pow(QoH0, 1 - bh);
                QTT = QoH0;
                QTT2 = Math.pow((QMA - 0.5 * QTT) / (AAH / ((1 - bh) * Ut)), 1 / (1 - bh));
                while(Math.abs(QTT2 - QTT) > err_limit_q) {
                    QTT = QTT2;
                    QTT2 = Math.pow((QMA - 0.5 * QTT) / (AAH / ((1 - bh) * Ut)), 1 / (1 - bh));
                }
                Q[i] = QTT2;
            } else if(bh > 1.0) {
                if (QoH0 > 0.0) {
                    SoH0 = AAH / (1 - bh) * (Math.pow(QoH0, 1 - bh) - Math.pow(0.01, 1 - bh));
                } else {
                    SoH0 = 0.0;
                }
                QMA = IAH - 0.5 * QoH0 + SoH0 + (AAH / ((1 - bh) * Ut)) * Math.pow(0.01, 1 - bh);
                QTT = QoH0;
                QTT2 = Math.pow((QMA - 0.5 * QTT) / (AAH / ((1 - bh) * Ut)), 1 / (1 - bh));
                while(Math.abs(QTT2 - QTT) > err_limit_q) {
                    QTT = QTT2;
                    QTT2 = Math.pow((QMA - 0.5 * QTT) / (AAH / ((1 - bh) * Ut)), 1 / (1 - bh));
                }
                Q[i] = QTT2;
            }
        } else { // 退水
            Q[i] = QoH0 * Math.pow(1 + bh / AAH * Math.pow(QoH0, bh), (-1) / bh);
        }
    }

    private List<String> toStringList(double[] array0){
        String[] array = new String[array0.length];
        for(int i = 0; i < array.length; i++){
            array[i] = String.valueOf(array0[i]);
        }
        return Arrays.asList(array);
    }

    @Override
    public void init() {

    }

    @Override
    public PredictResult predict() {
        PredictResult result = new PredictResult();
        Map<String, double[]> map = HBHDmodel();

        //河道出口流量
        double[] QSim = map.get("QSim");

        result.setQSim(map.get("QSim"));


        //更新状态
        List<String> QX = toStringList(Q0);

        Map<String, List<String>> status = new HashMap<>();
        status.put("QX", QX);

        result.setNewStatus(status);

        return result;
    }
}
